import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns; sns.set()
import pickle
import arviz as az
import pymc3 as pm
from matplotlib.colors import to_rgb
import scipy.stats as stats
from IPython.display import display
import matplotlib as mpl
%load_ext autoreload
%autoreload 2
import plotting_lib
writeOut = True
outPathPlots = "../plots/statistical_model_three_factors_filter_weak/"
outPathData = "../derived_data/statistical_model_three_factors_filter_weak/"
widthMM = 190
widthInch = widthMM / 25.4
ratio = 0.66666
heigthInch = ratio*widthInch
SMALL_SIZE = 8
MEDIUM_SIZE = 10
BIGGER_SIZE = 12
plt.rc('font', size=SMALL_SIZE) # controls default text sizes
plt.rc('axes', titlesize=SMALL_SIZE) # fontsize of the axes title
plt.rc('axes', labelsize=MEDIUM_SIZE) # fontsize of the x and y labels
plt.rc('xtick', labelsize=SMALL_SIZE) # fontsize of the tick labels
plt.rc('ytick', labelsize=SMALL_SIZE) # fontsize of the tick labels
plt.rc('legend', fontsize=SMALL_SIZE) # legend fontsize
plt.rc('figure', titlesize=BIGGER_SIZE) # fontsize of the figure title
sns.set_style("ticks")
dpi = 300
sizes = [SMALL_SIZE,MEDIUM_SIZE,BIGGER_SIZE]
numSamples = 1000
numCores = 4
numTune = 1000
numPredSamples = 2000
random_seed=36514535
target_accept = 0.99
datafile = "../derived_data/preprocessing/preprocessed_filter_weak.dat"
with open(datafile, "rb") as f:
x1,x2,x3,df,dataZ,dictMeanStd,dictTreatment,dictSoftware = pickle.load(f)
Show that everything is correct:
display(pd.DataFrame.from_dict({'x1':x1,'x2':x2,'x3':x3}))
| x1 | x2 | x3 | |
|---|---|---|---|
| 0 | 0 | 5 | 115 |
| 1 | 1 | 5 | 115 |
| 2 | 0 | 5 | 116 |
| 3 | 1 | 5 | 116 |
| 4 | 0 | 5 | 117 |
| ... | ... | ... | ... |
| 273 | 1 | 9 | 52 |
| 274 | 0 | 9 | 53 |
| 275 | 1 | 9 | 53 |
| 276 | 0 | 9 | 54 |
| 277 | 1 | 9 | 54 |
278 rows × 3 columns
x1 indicates the software used, x2 indicates the treatment applied and x3 the sample.
for surfaceParam,(mean,std) in dictMeanStd.items():
print("Surface parameter {} has mean {} and standard deviation {}".format(surfaceParam,mean,std))
Surface parameter epLsar has mean 0.0032388209205305753 and standard deviation 0.0019378273835719989 Surface parameter Rsquared has mean 0.9974096825435252 and standard deviation 0.007283582118542012 Surface parameter Asfc has mean 14.919474245449283 and standard deviation 12.47068676838922 Surface parameter Smfc has mean 1.155270960424856 and standard deviation 7.13503174525663 Surface parameter HAsfc9 has mean 0.44593694325514915 and standard deviation 0.7912033512620836 Surface parameter HAsfc81 has mean 0.9300206156734742 and standard deviation 2.3638534390774013
for k,v in sorted(dictTreatment.items(), key=lambda x: x[0]):
print("Number {} encodes treatment {}".format(k,v))
Number 0 encodes treatment BrushDirt Number 1 encodes treatment BrushNoDirt Number 2 encodes treatment Clover Number 3 encodes treatment Clover+dust Number 4 encodes treatment Control Number 5 encodes treatment Dry bamboo Number 6 encodes treatment Dry grass Number 7 encodes treatment Dry lucerne Number 8 encodes treatment Grass Number 9 encodes treatment Grass+dust Number 10 encodes treatment RubDirt
for k,v in dictSoftware.items():
print("Number {} encodes software {}".format(k,v))
Number 0 encodes software ConfoMap Number 1 encodes software Toothfrax
display(dataZ)
| index | TreatmentNumber | SoftwareNumber | DatasetNumber | NameNumber | epLsar_z | Rsquared_z | Asfc_z | Smfc_z | HAsfc9_z | HAsfc81_z | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 0 | 5 | 0 | 0 | 115 | 0.608031 | 0.081494 | -0.261684 | -0.120632 | -0.391977 | -0.239736 |
| 1 | 1 | 5 | 1 | 0 | 115 | 0.764866 | 0.295228 | -0.368764 | -0.145206 | -0.392397 | -0.240365 |
| 2 | 2 | 5 | 0 | 0 | 116 | 1.355641 | -0.166422 | 0.043912 | -0.120632 | -0.346351 | -0.268091 |
| 3 | 3 | 5 | 1 | 0 | 116 | 1.350574 | 0.282460 | -0.137943 | -0.145206 | -0.349727 | -0.282929 |
| 4 | 4 | 5 | 0 | 0 | 117 | 0.930308 | -0.359987 | -0.137793 | -0.120632 | -0.233444 | -0.221925 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 273 | 273 | 9 | 1 | 2 | 52 | 0.611602 | 0.267769 | -1.050437 | 0.035516 | 0.493708 | 0.076860 |
| 274 | 274 | 9 | 0 | 2 | 53 | 0.084569 | 0.197735 | -0.966638 | -0.093723 | 0.242115 | 0.257597 |
| 275 | 275 | 9 | 1 | 2 | 53 | -0.051512 | 0.319804 | -0.975181 | -0.143224 | 0.644288 | 0.381453 |
| 276 | 276 | 9 | 0 | 2 | 54 | -1.041990 | 0.284041 | -1.077552 | 0.011489 | -0.095103 | -0.053253 |
| 277 | 277 | 9 | 1 | 2 | 54 | -1.308590 | 0.336005 | -1.081522 | -0.104672 | -0.126169 | -0.038105 |
278 rows × 11 columns
display(df)
| Dataset | Name | Software | Diet | Treatment | Before.after | NMP | NMP_cat | epLsar | Rsquared | Asfc | Smfc | HAsfc9 | HAsfc81 | NewEplsar | TreatmentNumber | SoftwareNumber | DatasetNumber | NameNumber | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | GuineaPigs | capor_2CC6B1_txP4_#1_1_100xL_1 | ConfoMap | Dry bamboo | Dry bamboo | NaN | 0.717312 | 0-5% | 0.004417 | 0.998003 | 11.656095 | 0.294557 | 0.135803 | 0.363319 | 0.019460 | 5 | 0 | 0 | 115 |
| 1 | GuineaPigs | capor_2CC6B1_txP4_#1_1_100xL_1 | Toothfrax | Dry bamboo | Dry bamboo | NaN | 0.717312 | 0-5% | 0.004721 | 0.999560 | 10.320730 | 0.119219 | 0.135471 | 0.361833 | NaN | 5 | 1 | 0 | 115 |
| 2 | GuineaPigs | capor_2CC6B1_txP4_#1_1_100xL_2 | ConfoMap | Dry bamboo | Dry bamboo | NaN | 1.674215 | 0-5% | 0.005866 | 0.996198 | 15.467083 | 0.294557 | 0.171903 | 0.296292 | 0.020079 | 5 | 0 | 0 | 116 |
| 3 | GuineaPigs | capor_2CC6B1_txP4_#1_1_100xL_2 | Toothfrax | Dry bamboo | Dry bamboo | NaN | 1.674215 | 0-5% | 0.005856 | 0.999467 | 13.199232 | 0.119219 | 0.169232 | 0.261217 | NaN | 5 | 1 | 0 | 116 |
| 4 | GuineaPigs | capor_2CC6B1_txP4_#1_1_100xL_3 | ConfoMap | Dry bamboo | Dry bamboo | NaN | 1.760409 | 0-5% | 0.005042 | 0.994788 | 13.201101 | 0.294557 | 0.261235 | 0.405422 | 0.019722 | 5 | 0 | 0 | 117 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 273 | Sheeps | L8-Ovis-90730-lm2sin-a | Toothfrax | Grass+dust | Grass+dust | NaN | 0.000000 | 0-5% | 0.004424 | 0.999360 | 1.819802 | 1.408678 | 0.836560 | 1.111706 | NaN | 9 | 1 | 2 | 52 |
| 274 | Sheeps | L8-Ovis-90764-lm2sin-a | ConfoMap | Grass+dust | Grass+dust | NaN | 0.000000 | 0-5% | 0.003403 | 0.998850 | 2.864831 | 0.486556 | 0.637499 | 1.538943 | 0.018978 | 9 | 0 | 2 | 53 |
| 275 | Sheeps | L8-Ovis-90764-lm2sin-a | Toothfrax | Grass+dust | Grass+dust | NaN | 0.000000 | 0-5% | 0.003139 | 0.999739 | 2.758297 | 0.133366 | 0.955699 | 1.831721 | NaN | 9 | 1 | 2 | 53 |
| 276 | Sheeps | L8-Ovis-90814-lm2sin-a | ConfoMap | Grass+dust | Grass+dust | NaN | 0.000000 | 0-5% | 0.001220 | 0.999479 | 1.481662 | 1.237247 | 0.370691 | 0.804138 | 0.017498 | 9 | 0 | 2 | 54 |
| 277 | Sheeps | L8-Ovis-90814-lm2sin-a | Toothfrax | Grass+dust | Grass+dust | NaN | 0.000000 | 0-5% | 0.000703 | 0.999857 | 1.432148 | 0.408433 | 0.346111 | 0.839946 | NaN | 9 | 1 | 2 | 54 |
278 rows × 19 columns
Prepare dicts for contrasts that are of interest later.
x1contrast_dict = {'{}_{}'.format(dictSoftware[0],dictSoftware[1]):[1,-1]}
x1contrast_dict
{'ConfoMap_Toothfrax': [1, -1]}
dfRawDiff = dataZ.groupby(['DatasetNumber','TreatmentNumber','NameNumber']).agg(list).applymap(lambda l : l[0]-l[1]).reset_index()
dfRawDiff
| DatasetNumber | TreatmentNumber | NameNumber | index | SoftwareNumber | epLsar_z | Rsquared_z | Asfc_z | Smfc_z | HAsfc9_z | HAsfc81_z | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 0 | 5 | 115 | -1 | -1 | -0.156835 | -0.213734 | 0.107080 | 0.024574 | 0.000420 | 0.000629 |
| 1 | 0 | 5 | 116 | -1 | -1 | 0.005067 | -0.448881 | 0.181855 | 0.024574 | 0.003376 | 0.014838 |
| 2 | 0 | 5 | 117 | -1 | -1 | 0.137060 | -0.615262 | 0.171673 | 0.024574 | 0.055406 | 0.020671 |
| 3 | 0 | 5 | 118 | -1 | -1 | 0.363489 | -0.380505 | 0.206489 | 0.029593 | 0.018505 | -0.035897 |
| 4 | 0 | 5 | 119 | -1 | -1 | 0.135927 | -0.441161 | 0.431840 | 0.029593 | 0.051468 | 0.035179 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 134 | 2 | 9 | 50 | -1 | -1 | 0.565233 | -0.166711 | 0.141982 | 0.016303 | 0.055662 | 0.046721 |
| 135 | 2 | 9 | 51 | -1 | -1 | 1.694319 | -0.156497 | 0.021729 | 0.016303 | 0.446416 | 0.141461 |
| 136 | 2 | 9 | 52 | -1 | -1 | -2.109816 | -0.348137 | 0.000667 | 0.390193 | -0.010465 | 0.004987 |
| 137 | 2 | 9 | 53 | -1 | -1 | 0.136080 | -0.122069 | 0.008543 | 0.049501 | -0.402172 | -0.123856 |
| 138 | 2 | 9 | 54 | -1 | -1 | 0.266600 | -0.051963 | 0.003970 | 0.116161 | 0.031066 | -0.015148 |
139 rows × 11 columns
variablesList = ['epLsar','Rsquared','Asfc','Smfc','HAsfc9']
for var in variablesList:
fig, axes = plt.subplots(1, 2,figsize=(12,6))
fig.suptitle('{}'.format(var))
sns.stripplot(data=dataZ,x='NameNumber',y='{}_z'.format(var),hue='SoftwareNumber',ax=axes[0]);
axes[0].set_title('Raw data')
sns.histplot(data=dfRawDiff,x='{}_z'.format(var),ax=axes[1]);
axes[1].set_title('Aggegrated difference by software'.format(var))
plt.show()
class ThreeFactorModel(pm.Model):
"""
Compute params of priors and hyperpriors.
"""
def getParams(self,x1,x2,x3,y):
# get lengths
Nx1Lvl = np.unique(x1).size
Nx2Lvl = np.unique(x2).size
Nx3Lvl = np.unique(x3).size
dims = (Nx1Lvl, Nx2Lvl, Nx3Lvl)
### get standard deviations
# convert to pandas dataframe to use their logic
df = pd.DataFrame.from_dict({'x1':x1,'x2':x2,'x3':x3,'y':y})
s1 = df.groupby('x1').std()['y'].max()
s2 = df.groupby('x2').std()['y'].max()
s3 = df.groupby('x3').std()['y'].max()
stdSingle = (s1, s2, s3)
return (dims, stdSingle)
def printParams(self,x1,x2,x3,y):
dims, stdSingle = self.getParams(x1,x2,x3,y)
Nx1Lvl, Nx2Lvl, Nx3Lvl = dims
s1, s2, s3 = stdSingle
print("The number of levels of the x variables are {}".format(dims))
print("The standard deviations used for the beta priors are {}".format(stdSingle))
def __init__(self,name,x1,x2,x3,y,model=None):
# call super's init first, passing model and name
super().__init__(name, model)
# get parameter of hyperpriors
dims, stdSingle = self.getParams(x1,x2,x3,y)
Nx1Lvl, Nx2Lvl, Nx3Lvl = dims
s1, s2, s3 = stdSingle
### hyperpriors ###
# observation hyperpriors
lamY = 1/30.
muGamma = 0.5
sigmaGamma = 2.
# prediction hyperpriors
sigma0 = pm.HalfNormal('sigma0',sd=1)
sigma1 = pm.HalfNormal('sigma1',sd=s1, shape=Nx1Lvl)
sigma2 = pm.HalfNormal('sigma2',sd=s2, shape=Nx2Lvl)
sigma3 = pm.HalfNormal('sigma3',sd=s3, shape=Nx3Lvl)
mu_b0 = pm.Normal('mu_b0', mu=0., sd=1)
mu_b1 = pm.Normal('mu_b1', mu=0., sd=1, shape=Nx1Lvl)
mu_b2 = pm.Normal('mu_b2', mu=0., sd=1, shape=Nx2Lvl)
mu_b3 = pm.Normal('mu_b3', mu=0., sd=1, shape=Nx3Lvl)
### priors ###
# observation priors
nuY = pm.Exponential('nuY',lam=lamY)
sigmaY = pm.Gamma('sigmaY',mu=muGamma, sigma=sigmaGamma)
# prediction priors
b0_dist = pm.Normal('b0_dist', mu=0, sd=1)
b0 = pm.Deterministic("b0", mu_b0 + b0_dist * sigma0)
b1_dist = pm.Normal('b1_dist', mu=0, sd=1)
b1 = pm.Deterministic("b1", mu_b1 + b1_dist * sigma1)
b2_dist = pm.Normal('b2_dist', mu=0, sd=1)
b2 = pm.Deterministic("b2", mu_b2 + b2_dist * sigma2)
b3_dist = pm.Normal('b3_dist', mu=0, sd=1)
b3 = pm.Deterministic("b3", mu_b3 + b3_dist * sigma3)
#### prediction ###
mu = pm.Deterministic('mu',b0 + b1[x1] + b2[x2] + b3[x3])
### observation ###
y = pm.StudentT('y',nu = nuY, mu=mu, sd=sigmaY, observed=y)
with pm.Model() as model:
epLsarModel = ThreeFactorModel('epLsar',x1,x2,x3,dataZ.epLsar_z.values)
epLsarModel.printParams(x1,x2,x3,dataZ.epLsar_z.values)
The number of levels of the x variables are (2, 11, 139) The standard deviations used for the beta priors are (1.0181033807478441, 1.4891257020651163, 1.4918655146500768)
try:
graph_epLsar = pm.model_to_graphviz(epLsarModel)
except:
graph_epLsar = "Could not make graph"
graph_epLsar
with epLsarModel as model:
prior_pred_epLsar = pm.sample_prior_predictive(samples=numPredSamples,random_seed=random_seed)
plotting_lib.plotPriorPredictive(widthInch,heigthInch,dpi,writeOut,outPathPlots,df,dictMeanStd,prior_pred_epLsar,dataZ.epLsar_z.values,'epLsar')
Prior choice is as intended: Broad over the data range.
with epLsarModel as model:
trace_epLsar = pm.sample(numSamples,cores=numCores,tune=numTune,max_treedepth=20, init='auto',target_accept=target_accept,random_seed=random_seed)
Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Multiprocess sampling (4 chains in 4 jobs) NUTS: [epLsar_b3_dist, epLsar_b2_dist, epLsar_b1_dist, epLsar_b0_dist, epLsar_sigmaY, epLsar_nuY, epLsar_mu_b3, epLsar_mu_b2, epLsar_mu_b1, epLsar_mu_b0, epLsar_sigma3, epLsar_sigma2, epLsar_sigma1, epLsar_sigma0]
Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 590 seconds. The estimated number of effective samples is smaller than 200 for some parameters.
with epLsarModel as model:
if writeOut:
with open(outPathData + 'model_{}.pkl'.format('epLsar'), 'wb') as buff:
pickle.dump({'model':epLsarModel, 'trace': trace_epLsar}, buff)
with epLsarModel as model:
dataTrace_epLsar = az.from_pymc3(trace=trace_epLsar)
pm.summary(dataTrace_epLsar,hdi_prob=0.95).round(2)
| mean | sd | hdi_2.5% | hdi_97.5% | mcse_mean | mcse_sd | ess_mean | ess_sd | ess_bulk | ess_tail | r_hat | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| epLsar_mu_b0 | -0.04 | 0.79 | -1.50 | 1.64 | 0.01 | 0.01 | 4578.0 | 2208.0 | 4588.0 | 3218.0 | 1.00 |
| epLsar_mu_b1[0] | 0.02 | 0.68 | -1.42 | 1.25 | 0.01 | 0.01 | 3638.0 | 2339.0 | 3636.0 | 3182.0 | 1.00 |
| epLsar_mu_b1[1] | -0.04 | 0.68 | -1.41 | 1.22 | 0.01 | 0.01 | 3590.0 | 2011.0 | 3596.0 | 2887.0 | 1.00 |
| epLsar_mu_b2[0] | -0.14 | 0.53 | -1.16 | 0.89 | 0.01 | 0.01 | 1993.0 | 1993.0 | 1996.0 | 2449.0 | 1.00 |
| epLsar_mu_b2[1] | -0.18 | 0.52 | -1.22 | 0.81 | 0.01 | 0.01 | 1963.0 | 1825.0 | 1958.0 | 2574.0 | 1.00 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| epLsar_mu[273] | 0.18 | 0.78 | -1.65 | 0.90 | 0.05 | 0.03 | 283.0 | 283.0 | 555.0 | 965.0 | 1.01 |
| epLsar_mu[274] | 0.05 | 0.11 | -0.17 | 0.26 | 0.00 | 0.00 | 4555.0 | 2784.0 | 4392.0 | 3415.0 | 1.00 |
| epLsar_mu[275] | -0.01 | 0.11 | -0.22 | 0.22 | 0.00 | 0.00 | 4624.0 | 2478.0 | 4478.0 | 3058.0 | 1.00 |
| epLsar_mu[276] | -1.11 | 0.14 | -1.38 | -0.87 | 0.00 | 0.00 | 4435.0 | 4420.0 | 4210.0 | 3317.0 | 1.00 |
| epLsar_mu[277] | -1.18 | 0.14 | -1.43 | -0.92 | 0.00 | 0.00 | 4418.0 | 4415.0 | 4205.0 | 3404.0 | 1.00 |
743 rows × 11 columns
plotting_lib.plotDiagnostics(widthInch,heigthInch,dpi,writeOut,outPathPlots,trace_epLsar,dataTrace_epLsar,'epLsar')
/home/bob/Documents/Projekt_Neuwied/SSFA/ssfa-env/lib/python3.7/site-packages/arviz/plots/backends/matplotlib/pairplot.py:216: UserWarning: rcParams['plot.max_subplots'] (40) is smaller than the number of resulting pair plots with these variables, generating only a 8x8 grid UserWarning,
with epLsarModel as model:
plotting_lib.plotTracesB(widthInch,heigthInch,dpi,writeOut,outPathPlots,trace_epLsar,'epLsar')
with epLsarModel as model:
plotting_lib.pm.energyplot(trace_epLsar)
with epLsarModel as model:
posterior_pred_epLsar = pm.sample_posterior_predictive(trace_epLsar,samples=numPredSamples,random_seed=random_seed)
/home/bob/Documents/Projekt_Neuwied/SSFA/ssfa-env/lib/python3.7/site-packages/pymc3/sampling.py:1708: UserWarning: samples parameter is smaller than nchains times ndraws, some draws and/or chains may not be represented in the returned posterior predictive sample "samples parameter is smaller than nchains times ndraws, some draws "
plotting_lib.plotPriorPosteriorPredictive(widthInch,heigthInch,dpi,writeOut,outPathPlots,df,dictMeanStd,prior_pred_epLsar,posterior_pred_epLsar,dataZ.epLsar_z.values,'epLsar')
with epLsarModel as model:
pm_data_epLsar = az.from_pymc3(trace=trace_epLsar,prior=prior_pred_epLsar,posterior_predictive=posterior_pred_epLsar)
arviz.data.io_pymc3 - WARNING - posterior predictive variable epLsar_y's shape not compatible with number of chains and draws. This can mean that some draws or even whole chains are not represented.
plotting_lib.plotPriorPosteriorB(widthInch,heigthInch,dpi,sizes,writeOut,outPathPlots,dictMeanStd,pm_data_epLsar,'epLsar')
plotting_lib.plotPosterior(widthInch,heigthInch,dpi,writeOut,outPathPlots,dictMeanStd,pm_data_epLsar,'epLsar')
plotting_lib.plotContrast(widthInch,heigthInch,dpi,writeOut,outPathPlots,dictMeanStd,x1contrast_dict,trace_epLsar,'epLsar')
with pm.Model() as model:
RsquaredModel = ThreeFactorModel('Rsquared',x1,x2,x3,dataZ["Rsquared_z"].values)
RsquaredModel.printParams(x1,x2,x3,dataZ["Rsquared_z"].values)
The number of levels of the x variables are (2, 11, 139) The standard deviations used for the beta priors are (1.364226512794052, 3.833663350609328, 8.30397573287845)
pm.model_to_graphviz(RsquaredModel)
with RsquaredModel as model:
prior_pred_Rsquared = pm.sample_prior_predictive(samples=numPredSamples,random_seed=random_seed)
plotting_lib.plotPriorPredictive(widthInch,heigthInch,dpi,writeOut,outPathPlots,df,dictMeanStd,prior_pred_Rsquared,dataZ["Rsquared_z"].values,'Rsquared')
with RsquaredModel as model:
trace_Rsquared = pm.sample(numSamples,cores=numCores,tune=numTune,max_treedepth=20, init='auto',target_accept=target_accept,random_seed=random_seed)
#fit_Rsquared = pm.fit(random_seed=random_seed)
#trace_Rsquared = fit_Rsquared.sample(draws=numSamples)
Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Multiprocess sampling (4 chains in 4 jobs) NUTS: [Rsquared_b3_dist, Rsquared_b2_dist, Rsquared_b1_dist, Rsquared_b0_dist, Rsquared_sigmaY, Rsquared_nuY, Rsquared_mu_b3, Rsquared_mu_b2, Rsquared_mu_b1, Rsquared_mu_b0, Rsquared_sigma3, Rsquared_sigma2, Rsquared_sigma1, Rsquared_sigma0]
Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 827 seconds. The estimated number of effective samples is smaller than 200 for some parameters.
with RsquaredModel as model:
if writeOut:
with open(outPathData + 'model_{}.pkl'.format('Rsquared'), 'wb') as buff:
pickle.dump({'model': RsquaredModel, 'trace': trace_Rsquared}, buff)
with RsquaredModel as model:
dataTrace_Rsquared = az.from_pymc3(trace=trace_Rsquared)
pm.summary(dataTrace_Rsquared,hdi_prob=0.95).round(2)
| mean | sd | hdi_2.5% | hdi_97.5% | mcse_mean | mcse_sd | ess_mean | ess_sd | ess_bulk | ess_tail | r_hat | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| Rsquared_mu_b0 | 0.05 | 0.82 | -1.58 | 1.64 | 0.01 | 0.01 | 3512.0 | 2594.0 | 3532.0 | 3278.0 | 1.0 |
| Rsquared_mu_b1[0] | -0.08 | 0.73 | -1.44 | 1.38 | 0.01 | 0.01 | 2577.0 | 2577.0 | 2580.0 | 2991.0 | 1.0 |
| Rsquared_mu_b1[1] | 0.11 | 0.73 | -1.31 | 1.57 | 0.01 | 0.01 | 2693.0 | 2356.0 | 2682.0 | 2826.0 | 1.0 |
| Rsquared_mu_b2[0] | 0.08 | 0.53 | -0.98 | 1.16 | 0.01 | 0.01 | 1410.0 | 1410.0 | 1396.0 | 1961.0 | 1.0 |
| Rsquared_mu_b2[1] | -0.11 | 0.51 | -1.04 | 0.98 | 0.02 | 0.01 | 1221.0 | 1221.0 | 1228.0 | 1654.0 | 1.0 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| Rsquared_mu[273] | 0.21 | 0.09 | 0.05 | 0.40 | 0.00 | 0.00 | 4436.0 | 3566.0 | 4200.0 | 3064.0 | 1.0 |
| Rsquared_mu[274] | 0.14 | 0.09 | -0.05 | 0.31 | 0.00 | 0.00 | 4024.0 | 2720.0 | 4023.0 | 2677.0 | 1.0 |
| Rsquared_mu[275] | 0.38 | 0.09 | 0.20 | 0.56 | 0.00 | 0.00 | 4009.0 | 3705.0 | 3956.0 | 2511.0 | 1.0 |
| Rsquared_mu[276] | 0.19 | 0.11 | -0.01 | 0.41 | 0.00 | 0.00 | 4329.0 | 2580.0 | 4159.0 | 3037.0 | 1.0 |
| Rsquared_mu[277] | 0.43 | 0.11 | 0.22 | 0.64 | 0.00 | 0.00 | 4296.0 | 3506.0 | 4111.0 | 2970.0 | 1.0 |
743 rows × 11 columns
plotting_lib.plotDiagnostics(widthInch,heigthInch,dpi,writeOut,outPathPlots,dataTrace_Rsquared,trace_Rsquared,'Rsquared')
/home/bob/Documents/Projekt_Neuwied/SSFA/ssfa-env/lib/python3.7/site-packages/arviz/data/io_pymc3.py:89: FutureWarning: Using `from_pymc3` without the model will be deprecated in a future release. Not using the model will return less accurate and less useful results. Make sure you use the model argument or call from_pymc3 within a model context. FutureWarning, /home/bob/Documents/Projekt_Neuwied/SSFA/ssfa-env/lib/python3.7/site-packages/arviz/data/io_pymc3.py:89: FutureWarning: Using `from_pymc3` without the model will be deprecated in a future release. Not using the model will return less accurate and less useful results. Make sure you use the model argument or call from_pymc3 within a model context. FutureWarning, /home/bob/Documents/Projekt_Neuwied/SSFA/ssfa-env/lib/python3.7/site-packages/arviz/plots/backends/matplotlib/pairplot.py:216: UserWarning: rcParams['plot.max_subplots'] (40) is smaller than the number of resulting pair plots with these variables, generating only a 8x8 grid UserWarning, /home/bob/Documents/Projekt_Neuwied/SSFA/ssfa-env/lib/python3.7/site-packages/arviz/data/io_pymc3.py:89: FutureWarning: Using `from_pymc3` without the model will be deprecated in a future release. Not using the model will return less accurate and less useful results. Make sure you use the model argument or call from_pymc3 within a model context. FutureWarning, /home/bob/Documents/Projekt_Neuwied/SSFA/ssfa-env/lib/python3.7/site-packages/arviz/data/io_pymc3.py:89: FutureWarning: Using `from_pymc3` without the model will be deprecated in a future release. Not using the model will return less accurate and less useful results. Make sure you use the model argument or call from_pymc3 within a model context. FutureWarning,
with RsquaredModel as model:
plotting_lib.plotTracesB(widthInch,heigthInch,dpi,writeOut,outPathPlots,trace_Rsquared,'Rsquared')
with RsquaredModel as model:
plotting_lib.pm.energyplot(trace_Rsquared)
with RsquaredModel as model:
posterior_pred_Rsquared = pm.sample_posterior_predictive(trace_Rsquared,samples=numPredSamples,random_seed=random_seed)
/home/bob/Documents/Projekt_Neuwied/SSFA/ssfa-env/lib/python3.7/site-packages/pymc3/sampling.py:1708: UserWarning: samples parameter is smaller than nchains times ndraws, some draws and/or chains may not be represented in the returned posterior predictive sample "samples parameter is smaller than nchains times ndraws, some draws "
plotting_lib.plotPriorPosteriorPredictive(widthInch,heigthInch,dpi,writeOut,outPathPlots,df,dictMeanStd,prior_pred_Rsquared,posterior_pred_Rsquared,dataZ["Rsquared_z"].values,'Rsquared')
with RsquaredModel as model:
pm_data_Rsquared = az.from_pymc3(trace=trace_Rsquared,prior=prior_pred_Rsquared,posterior_predictive=posterior_pred_Rsquared)
arviz.data.io_pymc3 - WARNING - posterior predictive variable Rsquared_y's shape not compatible with number of chains and draws. This can mean that some draws or even whole chains are not represented.
plotting_lib.plotPriorPosteriorB(widthInch,heigthInch,dpi,sizes,writeOut,outPathPlots,dictMeanStd,pm_data_Rsquared,'Rsquared')
plotting_lib.plotPosterior(widthInch,heigthInch,dpi,writeOut,outPathPlots,dictMeanStd,pm_data_Rsquared,'Rsquared')
plotting_lib.plotContrast(widthInch,heigthInch,dpi,writeOut,outPathPlots,dictMeanStd,x1contrast_dict,trace_Rsquared,'Rsquared')
with pm.Model() as model:
AsfcModel = ThreeFactorModel('Asfc',x1,x2,x3,dataZ["Asfc_z"].values)
AsfcModel.printParams(x1,x2,x3,dataZ["Asfc_z"].values)
The number of levels of the x variables are (2, 11, 139) The standard deviations used for the beta priors are (1.0834504037573414, 1.7347208877727194, 0.46758648719539714)
pm.model_to_graphviz(AsfcModel)
with AsfcModel as model:
prior_pred_Asfc = pm.sample_prior_predictive(samples=numPredSamples,random_seed=random_seed)
plotting_lib.plotPriorPredictive(widthInch,heigthInch,dpi,writeOut,outPathPlots,df,dictMeanStd,prior_pred_Asfc,dataZ["Asfc_z"].values,'Asfc')
Prior choice is as intended: Broad over the data range.
with AsfcModel as model:
trace_Asfc = pm.sample(numSamples,cores=numCores,tune=numTune,max_treedepth=20, init='auto',target_accept=target_accept,random_seed=random_seed)
#fit_Asfc = pm.fit(random_seed=random_seed)
#trace_Asfc = fit_Asfc.sample(draws=numSamples)
Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Multiprocess sampling (4 chains in 4 jobs) NUTS: [Asfc_b3_dist, Asfc_b2_dist, Asfc_b1_dist, Asfc_b0_dist, Asfc_sigmaY, Asfc_nuY, Asfc_mu_b3, Asfc_mu_b2, Asfc_mu_b1, Asfc_mu_b0, Asfc_sigma3, Asfc_sigma2, Asfc_sigma1, Asfc_sigma0]
Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 668 seconds. The rhat statistic is larger than 1.05 for some parameters. This indicates slight problems during sampling. The estimated number of effective samples is smaller than 200 for some parameters.
with AsfcModel as model:
if writeOut:
with open(outPathData + 'model_{}.pkl'.format('Asfc'), 'wb') as buff:
pickle.dump({'model': AsfcModel, 'trace': trace_Asfc}, buff)
with AsfcModel as model:
dataTrace_Asfc = az.from_pymc3(trace=trace_Asfc)
pm.summary(dataTrace_Asfc,hdi_prob=0.95).round(2)
| mean | sd | hdi_2.5% | hdi_97.5% | mcse_mean | mcse_sd | ess_mean | ess_sd | ess_bulk | ess_tail | r_hat | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| Asfc_mu_b0 | 0.05 | 0.83 | -1.55 | 1.72 | 0.02 | 0.01 | 2909.0 | 2255.0 | 2905.0 | 2923.0 | 1.0 |
| Asfc_mu_b1[0] | 0.12 | 0.71 | -1.35 | 1.45 | 0.01 | 0.01 | 2599.0 | 2599.0 | 2600.0 | 3001.0 | 1.0 |
| Asfc_mu_b1[1] | -0.07 | 0.71 | -1.55 | 1.28 | 0.01 | 0.01 | 2755.0 | 2673.0 | 2754.0 | 2993.0 | 1.0 |
| Asfc_mu_b2[0] | 1.09 | 0.58 | 0.02 | 2.26 | 0.01 | 0.01 | 1841.0 | 1841.0 | 1844.0 | 2312.0 | 1.0 |
| Asfc_mu_b2[1] | 0.87 | 0.55 | -0.20 | 1.97 | 0.02 | 0.01 | 1376.0 | 1376.0 | 1380.0 | 1877.0 | 1.0 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| Asfc_mu[273] | -1.16 | 0.10 | -1.35 | -0.96 | 0.00 | 0.00 | 4235.0 | 4235.0 | 4229.0 | 3084.0 | 1.0 |
| Asfc_mu[274] | -0.86 | 0.10 | -1.06 | -0.68 | 0.00 | 0.00 | 4135.0 | 4017.0 | 4147.0 | 2803.0 | 1.0 |
| Asfc_mu[275] | -1.08 | 0.10 | -1.28 | -0.90 | 0.00 | 0.00 | 4093.0 | 4048.0 | 4115.0 | 2683.0 | 1.0 |
| Asfc_mu[276] | -0.97 | 0.10 | -1.16 | -0.78 | 0.00 | 0.00 | 4176.0 | 4161.0 | 4171.0 | 2895.0 | 1.0 |
| Asfc_mu[277] | -1.19 | 0.10 | -1.38 | -1.00 | 0.00 | 0.00 | 4128.0 | 4117.0 | 4124.0 | 2964.0 | 1.0 |
743 rows × 11 columns
plotting_lib.plotDiagnostics(widthInch,heigthInch,dpi,writeOut,outPathPlots,trace_Asfc,dataTrace_Asfc,'Asfc')
/home/bob/Documents/Projekt_Neuwied/SSFA/ssfa-env/lib/python3.7/site-packages/arviz/plots/backends/matplotlib/pairplot.py:216: UserWarning: rcParams['plot.max_subplots'] (40) is smaller than the number of resulting pair plots with these variables, generating only a 8x8 grid UserWarning,
with AsfcModel as model:
plotting_lib.plotTracesB(widthInch,heigthInch,dpi,writeOut,outPathPlots,trace_Asfc,'Asfc')
with AsfcModel as model:
plotting_lib.pm.energyplot(trace_Asfc)
with AsfcModel as model:
posterior_pred_Asfc = pm.sample_posterior_predictive(trace_Asfc,samples=numPredSamples,random_seed=random_seed)
/home/bob/Documents/Projekt_Neuwied/SSFA/ssfa-env/lib/python3.7/site-packages/pymc3/sampling.py:1708: UserWarning: samples parameter is smaller than nchains times ndraws, some draws and/or chains may not be represented in the returned posterior predictive sample "samples parameter is smaller than nchains times ndraws, some draws "
plotting_lib.plotPriorPosteriorPredictive(widthInch,heigthInch,dpi,writeOut,outPathPlots,df,dictMeanStd,prior_pred_Asfc,posterior_pred_Asfc,dataZ["Asfc_z"].values,'Asfc')
with AsfcModel as model:
pm_data_Asfc = az.from_pymc3(trace=trace_Asfc,prior=prior_pred_Asfc,posterior_predictive=posterior_pred_Asfc)
arviz.data.io_pymc3 - WARNING - posterior predictive variable Asfc_y's shape not compatible with number of chains and draws. This can mean that some draws or even whole chains are not represented.
plotting_lib.plotPriorPosteriorB(widthInch,heigthInch,dpi,sizes,writeOut,outPathPlots,dictMeanStd,pm_data_Asfc,'Asfc')
plotting_lib.plotPosterior(widthInch,heigthInch,dpi,writeOut,outPathPlots,dictMeanStd,pm_data_Asfc,'Asfc')
plotting_lib.plotContrast(widthInch,heigthInch,dpi,writeOut,outPathPlots,dictMeanStd,x1contrast_dict,trace_Asfc,'Asfc')
with pm.Model() as model:
SmfcModel = ThreeFactorModel('Smfc',x1,x2,x3,dataZ.Smfc_z.values)
SmfcModel.printParams(x1,x2,x3,dataZ.Smfc_z.values)
The number of levels of the x variables are (2, 11, 139) The standard deviations used for the beta priors are (1.2903046708375803, 3.554940942907768, 10.68509404141528)
pm.model_to_graphviz(SmfcModel)
with SmfcModel as model:
prior_pred_Smfc = pm.sample_prior_predictive(samples=numPredSamples,random_seed=random_seed)
plotting_lib.plotPriorPredictive(widthInch,heigthInch,dpi,writeOut,outPathPlots,df,dictMeanStd,prior_pred_Smfc,dataZ.Smfc_z.values,'Smfc')
Prior choice is as intended: Broad over the data range.
with SmfcModel as model:
trace_Smfc = pm.sample(numSamples,cores=numCores,tune=numTune,max_treedepth=20, init='auto',target_accept=0.8,random_seed=random_seed)
#fit_Smfc = pm.fit(random_seed=random_seed)
#trace_Smfc = fit_Smfc.sample(draws=numSamples)
Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Multiprocess sampling (4 chains in 4 jobs) NUTS: [Smfc_b3_dist, Smfc_b2_dist, Smfc_b1_dist, Smfc_b0_dist, Smfc_sigmaY, Smfc_nuY, Smfc_mu_b3, Smfc_mu_b2, Smfc_mu_b1, Smfc_mu_b0, Smfc_sigma3, Smfc_sigma2, Smfc_sigma1, Smfc_sigma0]
--------------------------------------------------------------------------- KeyboardInterrupt Traceback (most recent call last) ~/Documents/Projekt_Neuwied/SSFA/ssfa-env/lib/python3.7/site-packages/pymc3/sampling.py in _mp_sample(draws, tune, step, chains, cores, chain, random_seed, start, progressbar, trace, model, callback, discard_tuned_samples, mp_ctx, pickle_backend, **kwargs) 1485 with sampler: -> 1486 for draw in sampler: 1487 trace = traces[draw.chain - chain] ~/Documents/Projekt_Neuwied/SSFA/ssfa-env/lib/python3.7/site-packages/pymc3/parallel_sampling.py in __iter__(self) 491 while self._active: --> 492 draw = ProcessAdapter.recv_draw(self._active) 493 proc, is_last, draw, tuning, stats, warns = draw ~/Documents/Projekt_Neuwied/SSFA/ssfa-env/lib/python3.7/site-packages/pymc3/parallel_sampling.py in recv_draw(processes, timeout) 351 pipes = [proc._msg_pipe for proc in processes] --> 352 ready = multiprocessing.connection.wait(pipes) 353 if not ready: /usr/lib/python3.7/multiprocessing/connection.py in wait(object_list, timeout) 919 while True: --> 920 ready = selector.select(timeout) 921 if ready: /usr/lib/python3.7/selectors.py in select(self, timeout) 414 try: --> 415 fd_event_list = self._selector.poll(timeout) 416 except InterruptedError: KeyboardInterrupt: During handling of the above exception, another exception occurred: ValueError Traceback (most recent call last) <ipython-input-85-755e53d074df> in <module> 1 with SmfcModel as model: ----> 2 trace_Smfc = pm.sample(numSamples,cores=numCores,tune=numTune,max_treedepth=20, init='auto',target_accept=0.8,random_seed=random_seed) 3 #fit_Smfc = pm.fit(random_seed=random_seed) 4 #trace_Smfc = fit_Smfc.sample(draws=numSamples) ~/Documents/Projekt_Neuwied/SSFA/ssfa-env/lib/python3.7/site-packages/pymc3/sampling.py in sample(draws, step, init, n_init, start, trace, chain_idx, chains, cores, tune, progressbar, model, random_seed, discard_tuned_samples, compute_convergence_checks, callback, return_inferencedata, idata_kwargs, mp_ctx, pickle_backend, **kwargs) 543 _print_step_hierarchy(step) 544 try: --> 545 trace = _mp_sample(**sample_args, **parallel_args) 546 except pickle.PickleError: 547 _log.warning("Could not pickle model, sampling singlethreaded.") ~/Documents/Projekt_Neuwied/SSFA/ssfa-env/lib/python3.7/site-packages/pymc3/sampling.py in _mp_sample(draws, tune, step, chains, cores, chain, random_seed, start, progressbar, trace, model, callback, discard_tuned_samples, mp_ctx, pickle_backend, **kwargs) 1510 except KeyboardInterrupt: 1511 if discard_tuned_samples: -> 1512 traces, length = _choose_chains(traces, tune) 1513 else: 1514 traces, length = _choose_chains(traces, 0) ~/Documents/Projekt_Neuwied/SSFA/ssfa-env/lib/python3.7/site-packages/pymc3/sampling.py in _choose_chains(traces, tune) 1528 lengths = [max(0, len(trace) - tune) for trace in traces] 1529 if not sum(lengths): -> 1530 raise ValueError("Not enough samples to build a trace.") 1531 1532 idxs = np.argsort(lengths)[::-1] ValueError: Not enough samples to build a trace.
Analysis stopped here because sampling did not converge. As the plot shows, some data points are very far away from the others, which would require the analysis to be based on more heavy-tailed distributions.
with pm.Model() as model:
HAsfc9Model = ThreeFactorModel('HAsfc9',x1,x2,x3,dataZ["HAsfc9_z"].values)
HAsfc9Model.printParams(x1,x2,x3,dataZ["HAsfc9_z"].values)
The number of levels of the x variables are (2, 11, 139) The standard deviations used for the beta priors are (1.2742199511468122, 3.5115369692183998, 5.975437122957168)
pm.model_to_graphviz(HAsfc9Model)
with HAsfc9Model as model:
prior_pred_HAsfc9 = pm.sample_prior_predictive(samples=numPredSamples,random_seed=random_seed)
plotting_lib.plotPriorPredictive(widthInch,heigthInch,dpi,writeOut,outPathPlots,df,dictMeanStd,prior_pred_HAsfc9,dataZ["HAsfc9_z"].values,'HAsfc9')
Prior choice is as intended: Broad over the data range.
with HAsfc9Model as model:
trace_HAsfc9 = pm.sample(numSamples,cores=numCores,tune=numTune,max_treedepth=20, init='auto',target_accept=target_accept,random_seed=random_seed)
#fit_HAsfc9 = pm.fit(random_seed=random_seed)
#trace_HAsfc9 = fit_HAsfc9.sample(draws=numSamples)
Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Multiprocess sampling (4 chains in 4 jobs) NUTS: [HAsfc9_b3_dist, HAsfc9_b2_dist, HAsfc9_b1_dist, HAsfc9_b0_dist, HAsfc9_sigmaY, HAsfc9_nuY, HAsfc9_mu_b3, HAsfc9_mu_b2, HAsfc9_mu_b1, HAsfc9_mu_b0, HAsfc9_sigma3, HAsfc9_sigma2, HAsfc9_sigma1, HAsfc9_sigma0]
Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 5014 seconds. The rhat statistic is larger than 1.05 for some parameters. This indicates slight problems during sampling. The estimated number of effective samples is smaller than 200 for some parameters.
with HAsfc9Model as model:
if writeOut:
with open(outPathData + 'model_{}.pkl'.format('HAsfc9'), 'wb') as buff:
pickle.dump({'model': HAsfc9Model, 'trace': trace_HAsfc9}, buff)
with HAsfc9Model as model:
dataTrace_HAsfc9 = az.from_pymc3(trace=trace_HAsfc9)
pm.summary(dataTrace_HAsfc9,hdi_prob=0.95).round(2)
| mean | sd | hdi_2.5% | hdi_97.5% | mcse_mean | mcse_sd | ess_mean | ess_sd | ess_bulk | ess_tail | r_hat | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| HAsfc9_mu_b0 | -0.02 | 0.83 | -1.68 | 1.58 | 0.01 | 0.01 | 4930.0 | 1763.0 | 4905.0 | 3032.0 | 1.00 |
| HAsfc9_mu_b1[0] | 0.01 | 0.71 | -1.38 | 1.45 | 0.01 | 0.01 | 4889.0 | 2383.0 | 4870.0 | 3394.0 | 1.00 |
| HAsfc9_mu_b1[1] | -0.01 | 0.71 | -1.34 | 1.51 | 0.01 | 0.01 | 4937.0 | 2152.0 | 4934.0 | 3077.0 | 1.00 |
| HAsfc9_mu_b2[0] | -0.26 | 0.53 | -1.37 | 0.74 | 0.01 | 0.01 | 1869.0 | 1869.0 | 1857.0 | 2582.0 | 1.00 |
| HAsfc9_mu_b2[1] | 0.07 | 0.52 | -0.99 | 1.02 | 0.01 | 0.01 | 1892.0 | 1892.0 | 1902.0 | 2461.0 | 1.00 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| HAsfc9_mu[273] | 0.47 | 0.03 | 0.42 | 0.52 | 0.00 | 0.00 | 3427.0 | 3302.0 | 3426.0 | 3592.0 | 1.00 |
| HAsfc9_mu[274] | 0.46 | 0.22 | 0.20 | 0.73 | 0.02 | 0.01 | 159.0 | 159.0 | 248.0 | 2818.0 | 1.01 |
| HAsfc9_mu[275] | 0.42 | 0.22 | 0.16 | 0.69 | 0.02 | 0.01 | 159.0 | 159.0 | 247.0 | 2869.0 | 1.01 |
| HAsfc9_mu[276] | -0.09 | 0.02 | -0.12 | -0.06 | 0.00 | 0.00 | 3060.0 | 2461.0 | 3936.0 | 3005.0 | 1.00 |
| HAsfc9_mu[277] | -0.13 | 0.02 | -0.16 | -0.10 | 0.00 | 0.00 | 3008.0 | 2507.0 | 3864.0 | 3304.0 | 1.00 |
743 rows × 11 columns
plotting_lib.plotDiagnostics(widthInch,heigthInch,dpi,writeOut,outPathPlots,trace_HAsfc9,dataTrace_HAsfc9,'HAsfc9')
/home/bob/Documents/Projekt_Neuwied/SSFA/ssfa-env/lib/python3.7/site-packages/arviz/plots/backends/matplotlib/pairplot.py:216: UserWarning: rcParams['plot.max_subplots'] (40) is smaller than the number of resulting pair plots with these variables, generating only a 8x8 grid UserWarning,
with HAsfc9Model as model:
plotting_lib.plotTracesB(widthInch,heigthInch,dpi,writeOut,outPathPlots,trace_HAsfc9,'HAsfc9')
with HAsfc9Model as model:
plotting_lib.pm.energyplot(trace_HAsfc9)
with HAsfc9Model as model:
posterior_pred_HAsfc9 = pm.sample_posterior_predictive(trace_HAsfc9,samples=numPredSamples,random_seed=random_seed)
/home/bob/Documents/Projekt_Neuwied/SSFA/ssfa-env/lib/python3.7/site-packages/pymc3/sampling.py:1708: UserWarning: samples parameter is smaller than nchains times ndraws, some draws and/or chains may not be represented in the returned posterior predictive sample "samples parameter is smaller than nchains times ndraws, some draws "
plotting_lib.plotPriorPosteriorPredictive(widthInch,heigthInch,dpi,writeOut,outPathPlots,df,dictMeanStd,prior_pred_HAsfc9,posterior_pred_HAsfc9,dataZ["HAsfc9_z"].values,'HAsfc9')
with HAsfc9Model as model:
pm_data_HAsfc9 = az.from_pymc3(trace=trace_HAsfc9,prior=prior_pred_HAsfc9,posterior_predictive=posterior_pred_HAsfc9)
arviz.data.io_pymc3 - WARNING - posterior predictive variable HAsfc9_y's shape not compatible with number of chains and draws. This can mean that some draws or even whole chains are not represented.
plotting_lib.plotPriorPosteriorB(widthInch,heigthInch,dpi,sizes,writeOut,outPathPlots,dictMeanStd,pm_data_HAsfc9,'HAsfc9')
plotting_lib.plotPosterior(widthInch,heigthInch,dpi,writeOut,outPathPlots,dictMeanStd,pm_data_HAsfc9,'HAsfc9')
plotting_lib.plotContrast(widthInch,heigthInch,dpi,writeOut,outPathPlots,dictMeanStd,x1contrast_dict,trace_HAsfc9,'HAsfc9')
The contrast plots between the two software packages are shown below again for each variable except Smfc
plotting_lib.plotContrast(widthInch,heigthInch,dpi,writeOut,outPathPlots,dictMeanStd,x1contrast_dict,trace_epLsar,'epLsar')
plotting_lib.plotContrast(widthInch,heigthInch,dpi,writeOut,outPathPlots,dictMeanStd,x1contrast_dict,trace_Rsquared,'Rsquared')
plotting_lib.plotContrast(widthInch,heigthInch,dpi,writeOut,outPathPlots,dictMeanStd,x1contrast_dict,trace_Asfc,'Asfc')
plotting_lib.plotContrast(widthInch,heigthInch,dpi,writeOut,outPathPlots,dictMeanStd,x1contrast_dict,trace_HAsfc9,'HAsfc9')
!jupyter nbconvert --to html Statistical_Model_ThreeFactor_filter_weak.ipynb
[NbConvertApp] Converting notebook Statistical_Model_ThreeFactor_filter_weak.ipynb to html [NbConvertApp] Writing 7303955 bytes to Statistical_Model_ThreeFactor_filter_weak.html
!jupyter nbconvert --to markdown Statistical_Model_ThreeFactor_filter_weak.ipynb
[NbConvertApp] Converting notebook Statistical_Model_ThreeFactor_filter_weak.ipynb to markdown [NbConvertApp] Support files will be in Statistical_Model_ThreeFactor_filter_weak_files/ [NbConvertApp] Making directory Statistical_Model_ThreeFactor_filter_weak_files [NbConvertApp] Making directory Statistical_Model_ThreeFactor_filter_weak_files [NbConvertApp] Making directory Statistical_Model_ThreeFactor_filter_weak_files [NbConvertApp] Making directory Statistical_Model_ThreeFactor_filter_weak_files [NbConvertApp] Making directory Statistical_Model_ThreeFactor_filter_weak_files [NbConvertApp] Making directory Statistical_Model_ThreeFactor_filter_weak_files [NbConvertApp] Making directory Statistical_Model_ThreeFactor_filter_weak_files [NbConvertApp] Making directory Statistical_Model_ThreeFactor_filter_weak_files [NbConvertApp] Making directory Statistical_Model_ThreeFactor_filter_weak_files [NbConvertApp] Making directory Statistical_Model_ThreeFactor_filter_weak_files [NbConvertApp] Making directory Statistical_Model_ThreeFactor_filter_weak_files [NbConvertApp] Making directory Statistical_Model_ThreeFactor_filter_weak_files [NbConvertApp] Making directory Statistical_Model_ThreeFactor_filter_weak_files [NbConvertApp] Making directory Statistical_Model_ThreeFactor_filter_weak_files [NbConvertApp] Making directory Statistical_Model_ThreeFactor_filter_weak_files [NbConvertApp] Making directory Statistical_Model_ThreeFactor_filter_weak_files [NbConvertApp] Making directory Statistical_Model_ThreeFactor_filter_weak_files [NbConvertApp] Making directory Statistical_Model_ThreeFactor_filter_weak_files [NbConvertApp] Making directory Statistical_Model_ThreeFactor_filter_weak_files [NbConvertApp] Making directory Statistical_Model_ThreeFactor_filter_weak_files [NbConvertApp] Making directory Statistical_Model_ThreeFactor_filter_weak_files [NbConvertApp] Making directory Statistical_Model_ThreeFactor_filter_weak_files [NbConvertApp] Making directory Statistical_Model_ThreeFactor_filter_weak_files [NbConvertApp] Making directory Statistical_Model_ThreeFactor_filter_weak_files [NbConvertApp] Making directory Statistical_Model_ThreeFactor_filter_weak_files [NbConvertApp] Making directory Statistical_Model_ThreeFactor_filter_weak_files [NbConvertApp] Making directory Statistical_Model_ThreeFactor_filter_weak_files [NbConvertApp] Making directory Statistical_Model_ThreeFactor_filter_weak_files [NbConvertApp] Making directory Statistical_Model_ThreeFactor_filter_weak_files [NbConvertApp] Making directory Statistical_Model_ThreeFactor_filter_weak_files [NbConvertApp] Making directory Statistical_Model_ThreeFactor_filter_weak_files [NbConvertApp] Making directory Statistical_Model_ThreeFactor_filter_weak_files [NbConvertApp] Making directory Statistical_Model_ThreeFactor_filter_weak_files [NbConvertApp] Making directory Statistical_Model_ThreeFactor_filter_weak_files [NbConvertApp] Making directory Statistical_Model_ThreeFactor_filter_weak_files [NbConvertApp] Making directory Statistical_Model_ThreeFactor_filter_weak_files [NbConvertApp] Making directory Statistical_Model_ThreeFactor_filter_weak_files [NbConvertApp] Making directory Statistical_Model_ThreeFactor_filter_weak_files [NbConvertApp] Making directory Statistical_Model_ThreeFactor_filter_weak_files [NbConvertApp] Making directory Statistical_Model_ThreeFactor_filter_weak_files [NbConvertApp] Making directory Statistical_Model_ThreeFactor_filter_weak_files [NbConvertApp] Making directory Statistical_Model_ThreeFactor_filter_weak_files [NbConvertApp] Making directory Statistical_Model_ThreeFactor_filter_weak_files [NbConvertApp] Making directory Statistical_Model_ThreeFactor_filter_weak_files [NbConvertApp] Making directory Statistical_Model_ThreeFactor_filter_weak_files [NbConvertApp] Writing 69184 bytes to Statistical_Model_ThreeFactor_filter_weak.md